1 Operaciones espaciales

Hay una enorma cantidad de operaciones que se le pueden hacer a los objetos geográficos. Calcular el área, el perímetro, el largo, distancias entre objetos, intersecciones, uniones, diferencias, realizar transformaciones y más! Estas manipulaciones son sumamente útiles ya que nos permitirán comprender mejor el fenómeno que estamos analizando.

Para más información de todas las operaciones posibles (o las principales) sugerimos revisar el capítulo 4 del libro Geocomputation with R.

2 Área, perímetro e intersección

Como siempre comenzamos cargando las librerías que vamos a usar. Recuerden que si no las tienen instaladas aún lo pueden hacer con install.packages(“NOMBRE_LIBRERIA”).

De paso cargamos un dataframe geográfico (formato .shp) que tiene los polígonos de las provincias de Argentina.

library(sf) 
library(tidyverse)
library(units) 

options(scipen = 999)

prov <- read_sf("https://raw.githubusercontent.com/martoalalu/clase-geo-salud/clase-2021/data/ar_provincias.geojson")

Veamos qué información tiene.

head(prov)

## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -72.13788 ymin: -46.0069 xmax: -56.66736 ymax: -24.09314
## Geodetic CRS:  WGS 84
## # A tibble: 6 x 2
##   provincia                                                             geometry
##   <chr>                                                       <MULTIPOLYGON [°]>
## 1 Buenos Aires       (((-62.4578 -38.75891, -62.45403 -38.75514, -62.43847 -38.…
## 2 Catamarca          (((-64.82775 -29.56298, -64.9697 -29.58564, -64.90726 -29.…
## 3 Chaco              (((-60.9271 -27.99942, -61.71717 -27.99393, -61.71229 -25.…
## 4 Chubut             (((-66.36792 -45.04317, -66.37458 -45.05514, -66.36292 -45…
## 5 Ciudad de Buenos … (((-58.44724 -34.69322, -58.46166 -34.70905, -58.52893 -34…
## 6 Córdoba            (((-63.39265 -34.99844, -65.10672 -34.99768, -65.10981 -33…

Ok, sólo el nombre y el campo de geometría, obvio. Hagamos un mapa!

ggplot(prov) +
    geom_sf()

Podemos rellenar cada provincia con un color distinto.

ggplot(prov) +
    geom_sf(aes(fill=provincia))

Ok. Ahora sí, empecemos a agregarle información a este dataframe. Calcular el área es muy simple con la función st_area de la librería sf.

st_area(prov) #nos devuelve en unidades raras, pasamos a km2

## Units: [m^2]
##  [1] 306889486297 101887635888  99886544213 224258973324    209979717
##  [6] 164683302391  88973340510  77861018534  75596524598  53192452823
## [11] 142688382924  90990885271 148787352879  29982817814  94604326333
## [16] 202590178533 155322607325  88301731606  75855929854 243981330500
## [21] 132921440646 136871676316  20909206585  22537877676

Mmm suenan raros estos números. Pasemos la unidad de medición a kilómetros cuadrados!

set_units(st_area(prov), km^2)

## Units: [km^2]
##  [1] 306889.4863 101887.6359  99886.5442 224258.9733    209.9797 164683.3024
##  [7]  88973.3405  77861.0185  75596.5246  53192.4528 142688.3829  90990.8853
## [13] 148787.3529  29982.8178  94604.3263 202590.1785 155322.6073  88301.7316
## [19]  75855.9299 243981.3305 132921.4406 136871.6763  20909.2066  22537.8777

Ahí va mejor. Pero seguimos sin saber a qué provincia corresponde cada registro. Lo pasamos a formato numérico y agregamos como columna nueva!

prov$area <- as.numeric(set_units(st_area(prov), km^2))
prov

## Simple feature collection with 24 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.57778 ymin: -55.05736 xmax: -53.59184 ymax: -21.77855
## Geodetic CRS:  WGS 84
## # A tibble: 24 x 3
##    provincia                                                     geometry   area
##  * <chr>                                               <MULTIPOLYGON [°]>  <dbl>
##  1 Buenos Aires     (((-62.4578 -38.75891, -62.45403 -38.75514, -62.4384… 3.07e5
##  2 Catamarca        (((-64.82775 -29.56298, -64.9697 -29.58564, -64.9072… 1.02e5
##  3 Chaco            (((-60.9271 -27.99942, -61.71717 -27.99393, -61.7122… 9.99e4
##  4 Chubut           (((-66.36792 -45.04317, -66.37458 -45.05514, -66.362… 2.24e5
##  5 Ciudad de Bueno… (((-58.44724 -34.69322, -58.46166 -34.70905, -58.528… 2.10e2
##  6 Córdoba          (((-63.39265 -34.99844, -65.10672 -34.99768, -65.109… 1.65e5
##  7 Corrientes       (((-57.78307 -30.43265, -57.84768 -30.47401, -57.867… 8.90e4
##  8 Entre Ríos       (((-58.5247 -33.46152, -58.53116 -33.49211, -58.5254… 7.79e4
##  9 Formosa          (((-58.24602 -26.73134, -58.25027 -26.75323, -58.261… 7.56e4
## 10 Jujuy            (((-64.83212 -24.47771, -64.86935 -24.5232, -64.8812… 5.32e4
## # … with 14 more rows

Bien! Ahora que tenemos una columna con el área podemos hacer un mapa coloreando cada provincia según los km2.

ggplot(prov) +
  geom_sf(aes(fill=area))

Mmm está mal, pero no tan mal. Lo que ya sabíamos, la Provincia de Buenos Aires es la más grande del país, pero también se ve cierto patrón, a medida que pasamos de Norte a Sur el gradiente se ve más claro, o sea que las provincias suelen ser más grandes…

Hagamos un gráfico de barras para ver este patrón.

ggplot(prov) + 
  geom_bar(aes(x= reorder(provincia, area), weight=area, fill=area)) +
  coord_flip()

Claro, al pasar del mapa al gráfico de barras perdemos la referencia geográfica. Agreguemos una columna con la región de la provincia para ver si las de la Patagonia suelen ser más grandes que las del Norte.

prov <- prov %>% 
  mutate(region =
  case_when(
        provincia %in% (c("Buenos Aires", "Córdoba", "La Pampa", "Entre Ríos", "Santa Fe", "Ciudad de Buenos Aires")) ~ "Pampeana",
        provincia %in% (c("Mendoza", "San Luis", "San Juan")) ~ "Cuyo",
        provincia %in% (c("La Rioja", "Catamarca", "Jujuy", "Salta", "Tucumán", "Santiago del Estero")) ~ "Noroeste",
        provincia %in% (c("Corrientes", "Misiones", "Formosa", "Chaco")) ~ "Noreste",
        TRUE ~ "Patagonia"
        )
  )

ggplot(prov) + 
  geom_bar(aes(x= reorder(provincia, area), weight=area, fill=region)) +
  coord_flip()

Ahí va mejor! Hay cierto patrón, 3 de las 4 provincias más grande son de la Patagonia y 5 de las 6 más chicas del Norte! Esto nos sirve para entender que al hacer análisis espacial a veces salir del mapa puede ser sumamente útil para ver patrones.

Una más y listo. Agreguemos una línea vertical con la media del área de las provincias para ver cuáles están por encima y cuáles por debajo.

ggplot(prov) + 
  geom_bar(aes(x= reorder(provincia, area), weight=area, fill=region)) +
 geom_hline(yintercept=mean(prov$area), linetype="dashed", 
                color = "black", size=0.5) +
  coord_flip()

¿Qué conclusiones sacan ahora?

Bueno, sigamos con el perímetro. La función es st_length.

st_length(prov)

## Units: [m]
##  [1] 4502992.56 2045513.65 1787796.08 3069110.91   61337.54 1887257.68
##  [7] 1523652.60 1371928.96 1823365.32 1356250.58 1801249.30 1592326.16
## [13] 2012790.78 1218772.88 1767529.09 2570198.77 2869031.24 1716796.00
## [19] 1364124.16 2773527.11 2143460.03 1701827.16 1449042.72  806231.04

Vemos que la unidad de medida son metros, lo pasamos a km con solo dividirlo por 1000 y agregamos a columna nueva.

prov$perimetro <- as.numeric(st_length(prov) / 1000)
prov

## Simple feature collection with 24 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.57778 ymin: -55.05736 xmax: -53.59184 ymax: -21.77855
## Geodetic CRS:  WGS 84
## # A tibble: 24 x 5
##    provincia                                    geometry   area region perimetro
##  * <chr>                              <MULTIPOLYGON [°]>  <dbl> <chr>      <dbl>
##  1 Buenos Aires  (((-62.4578 -38.75891, -62.45403 -38.7… 3.07e5 Pampe…    4503. 
##  2 Catamarca     (((-64.82775 -29.56298, -64.9697 -29.5… 1.02e5 Noroe…    2046. 
##  3 Chaco         (((-60.9271 -27.99942, -61.71717 -27.9… 9.99e4 Nores…    1788. 
##  4 Chubut        (((-66.36792 -45.04317, -66.37458 -45.… 2.24e5 Patag…    3069. 
##  5 Ciudad de Bu… (((-58.44724 -34.69322, -58.46166 -34.… 2.10e2 Pampe…      61.3
##  6 Córdoba       (((-63.39265 -34.99844, -65.10672 -34.… 1.65e5 Pampe…    1887. 
##  7 Corrientes    (((-57.78307 -30.43265, -57.84768 -30.… 8.90e4 Nores…    1524. 
##  8 Entre Ríos    (((-58.5247 -33.46152, -58.53116 -33.4… 7.79e4 Pampe…    1372. 
##  9 Formosa       (((-58.24602 -26.73134, -58.25027 -26.… 7.56e4 Nores…    1823. 
## 10 Jujuy         (((-64.83212 -24.47771, -64.86935 -24.… 5.32e4 Noroe…    1356. 
## # … with 14 more rows

Excelente! Ahora sí, al mapa.

ggplot(prov) +
  geom_sf(aes(fill=perimetro))

Y al gráfico de barras.

ggplot(prov) + 
  geom_bar(aes(x= reorder(provincia, perimetro), weight=perimetro, fill=region)) +
 geom_hline(yintercept=mean(prov$perimetro), linetype="dashed", 
                color = "black", size=0.5) +
  coord_flip()

3 Intersecciones

Otra de las operaciones espaciales más útiles es la de calcular intersecciones entre 2 objetos espaciales. Por ejemplo ver qué puntos (escuelas) caen dentro de un polígono (ej. barrio). Muchas veces la información disponible es de un universo mucho mayor a nuestro fenómeno de estudio y necesitamos hacer un recorte, ahí es dónde st_intersects nos puede ayudar.

Asimismo por ahora estuvimos manejando un tipo de dato geográfico (polígonos) y haciendo un estilo de mapas (coropléticos), probemos usando datos que sean puntos y hagamos otros tipos de mapas.

Vamos a usar un dataset disponibilizado por el Ministerio de Salud del Gobierno Nacional que contiene la ubicación exacta todos los efectores de salud del país.

#Cargamos los datos
salud <- read.csv("https://github.com/martoalalu/clase-geo-salud/raw/master/data/refes-hospitales.csv", fileEncoding = 'UTF-8')

#Vemos la clase de objeto que es
class(salud)

## [1] "data.frame"

head(salud)

##           CODIGO GLN                                                 NOMBRE
## 1 50380842150137  NA                                    PUESTO EL TORO (17)
## 2 50380632150075  NA PUESTO BARRO NEGRO (9) (EX LOTE ROSARIO DE RIO GRANDE)
## 3 50260212329179  NA                                              KIN SPORT
## 4 51260352329290  NA                   LABORATORIO ANALISIS CLINICOS GEROSA
## 5 51260352329289  NA                          LABORATORIO ANALISIS CLINICOS
## 6 51260352329293  NA             LABORATORIO ANATOMIA PATOLOGICO (SALVAGNO)
##   CODIGO_TIPOLOGIA
## 1           ESSIDT
## 2           ESSIDT
## 3           ESSIDT
## 4            ESSID
## 5            ESSID
## 6            ESSID
##                                                               TIPOLOGIA
## 1 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 2 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 3 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 4               Establecimiento de salud sin internación de diagnóstico
## 5               Establecimiento de salud sin internación de diagnóstico
## 6               Establecimiento de salud sin internación de diagnóstico
##   CODIGO_CATEGORIA_TIPOLOGIA
## 1                     ESSIDT
## 2                     ESSIDT
## 3                     ESSDIT
## 4                      ESSID
## 5                      ESSID
## 6                      ESSID
##                                                     CATEGORIA_TIPOLOGIA
## 1   Sin atención médica en forma periódica (menor a 3 veces por semana)
## 2          Con atención médica general por lo menos 3 días de la semana
## 3 Con atención médica diaria y con especialidades y/o otras profesiones
## 4                                      Laboratorio de Análisis Clínicos
## 5                                      Laboratorio de Análisis Clínicos
## 6                                    Laboratorio de Anatomía patológica
##   CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## 1                 21  Provincial               Público               38
## 2                 21  Provincial               Público               38
## 3                 23     Privado               Privado               26
## 4                 23     Privado               Privado               26
## 5                 23     Privado               Privado               26
## 6                 23     Privado               Privado               26
##   PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## 1     Jujuy                  84      Susques      38084030000
## 2     Jujuy                  63    San Pedro      38063090000
## 3    Chubut                  21    Escalante      26021030018
## 4    Chubut                  35    Futaleufú      26035030000
## 5    Chubut                  35    Futaleufú      26035030000
## 6    Chubut                  35    Futaleufú      26035030000
##            LOCALIDAD   CP                             DOMICILIO          TE1
## 1            EL TORO 4411 AVENIDA LAVALLE ESQUINA 24 DE OCTUBRE             
## 2        LA MENDIETA 4522 CARLOS GARDEL S/N° - LOTE BARRO NEGRO             
## 3 COMODORO RIVADAVIA 3240                 Avenida Rivadavia 633 0297-4475310
## 4             ESQUEL 9200                        25 de Mayo 498   2945452050
## 5             ESQUEL 9200                         Ameghino 1033 02945-452155
## 6             ESQUEL 9200                          Belgrano 347             
##            TE2 TE3 TE4   LATITUD  LONGITUD   ESTADO_GEO
## 1                      -23.20853 -66.82709 SINCONFIRMAR
## 2                      -24.30814 -64.93175 SINCONFIRMAR
## 3                      -45.86213 -67.48355   CONFIRMADO
## 4 294515503809         -42.91505 -71.31971   CONFIRMADO
## 5                      -42.91577 -71.31803   CONFIRMADO
## 6                      -42.91788 -71.31992 SINCONFIRMAR

summary(salud)

##      CODIGO                    GLN                   NOMBRE         
##  Min.   :10020012215200   Min.   :7798100490010   Length:28366      
##  1st Qu.:50063572202800   1st Qu.:9990454430010   Class :character  
##  Median :50500072360600   Median :9991301000120   Mode  :character  
##  Mean   :43797418878000   Mean   :9868266292960                     
##  3rd Qu.:51065682305200   3rd Qu.:9992051100060                     
##  Max.   :53940142795200   Max.   :9992574600000                     
##                           NA's   :28241                             
##  CODIGO_TIPOLOGIA    TIPOLOGIA         CODIGO_CATEGORIA_TIPOLOGIA
##  Length:28366       Length:28366       Length:28366              
##  Class :character   Class :character   Class :character          
##  Mode  :character   Mode  :character   Mode  :character          
##                                                                  
##                                                                  
##                                                                  
##                                                                  
##  CATEGORIA_TIPOLOGIA CODIGO_DEPENDENCIA DEPENDENCIA       
##  Length:28366        Min.   :20.00      Length:28366      
##  Class :character    1st Qu.:22.00      Class :character  
##  Mode  :character    Median :23.00      Mode  :character  
##                      Mean   :22.52                        
##                      3rd Qu.:23.00                        
##                      Max.   :32.00                        
##                                                           
##  ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA  PROVINCIA         CODIGO_DEPARTAMENTO
##  Length:28366          Min.   : 2.00    Length:28366       Min.   :  1.0      
##  Class :character      1st Qu.: 6.00    Class :character   1st Qu.: 28.0      
##  Mode  :character      Median :30.00    Mode  :character   Median : 70.0      
##                        Mean   :38.81                       Mean   :157.5      
##                        3rd Qu.:66.00                       3rd Qu.:140.0      
##                        Max.   :94.00                       Max.   :882.0      
##                                                                               
##  DEPARTAMENTO       CODIGO_LOCALIDAD       LOCALIDAD              CP           
##  Length:28366       Min.   :    6028010   Length:28366       Length:28366      
##  Class :character   1st Qu.: 6588100000   Class :character   Class :character  
##  Mode  :character   Median :30035045000   Mode  :character   Mode  :character  
##                     Mean   :37490327544                                        
##                     3rd Qu.:66028050000                                        
##                     Max.   :94014020000                                        
##                                                                                
##   DOMICILIO             TE1                TE2                TE3           
##  Length:28366       Length:28366       Length:28366       Length:28366      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##      TE4               LATITUD          LONGITUD       ESTADO_GEO       
##  Length:28366       Min.   :-54.83   Min.   :-72.89   Length:28366      
##  Class :character   1st Qu.:-34.70   1st Qu.:-65.35   Class :character  
##  Mode  :character   Median :-32.91   Median :-61.93   Mode  :character  
##                     Mean   :-32.52   Mean   :-62.47                     
##                     3rd Qu.:-28.47   3rd Qu.:-58.67                     
##                     Max.   :-21.94   Max.   :-53.65                     
##                     NA's   :2240     NA's   :2240

Lo cargamos como un dataframe normal no-geográfico pero al explorarlo vemos que tiene 2 columnas que sí precisan su ubicación exacta en la tierra, LATITUD y LONGITUD. Tenemos que hacerle entender a R que se trata de un dataframe geográfico, es fácil!

Primero filtramos y sólo nos quedamos con las observaciones que no son nulas ya que no podemos ubicar en un mapa objetos con latitud y longitud nulas. Luego con st_as_sf convertimos el dataframe a un objeto geográfico indicandole dónde están los datos de las coordenadas y le indicamos el sistema de referencia de coordenadas.

salud <- salud %>% 
  filter(!is.na(LONGITUD), !is.na(LATITUD)) %>% 
  st_as_sf(coords = c("LONGITUD", "LATITUD"), crs = 4326)

class(salud)

## [1] "sf"         "data.frame"

head(salud)

## Simple feature collection with 6 features and 23 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -71.31992 ymin: -45.86213 xmax: -64.93175 ymax: -23.20853
## Geodetic CRS:  WGS 84
##           CODIGO GLN                                                 NOMBRE
## 1 50380842150137  NA                                    PUESTO EL TORO (17)
## 2 50380632150075  NA PUESTO BARRO NEGRO (9) (EX LOTE ROSARIO DE RIO GRANDE)
## 3 50260212329179  NA                                              KIN SPORT
## 4 51260352329290  NA                   LABORATORIO ANALISIS CLINICOS GEROSA
## 5 51260352329289  NA                          LABORATORIO ANALISIS CLINICOS
## 6 51260352329293  NA             LABORATORIO ANATOMIA PATOLOGICO (SALVAGNO)
##   CODIGO_TIPOLOGIA
## 1           ESSIDT
## 2           ESSIDT
## 3           ESSIDT
## 4            ESSID
## 5            ESSID
## 6            ESSID
##                                                               TIPOLOGIA
## 1 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 2 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 3 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 4               Establecimiento de salud sin internación de diagnóstico
## 5               Establecimiento de salud sin internación de diagnóstico
## 6               Establecimiento de salud sin internación de diagnóstico
##   CODIGO_CATEGORIA_TIPOLOGIA
## 1                     ESSIDT
## 2                     ESSIDT
## 3                     ESSDIT
## 4                      ESSID
## 5                      ESSID
## 6                      ESSID
##                                                     CATEGORIA_TIPOLOGIA
## 1   Sin atención médica en forma periódica (menor a 3 veces por semana)
## 2          Con atención médica general por lo menos 3 días de la semana
## 3 Con atención médica diaria y con especialidades y/o otras profesiones
## 4                                      Laboratorio de Análisis Clínicos
## 5                                      Laboratorio de Análisis Clínicos
## 6                                    Laboratorio de Anatomía patológica
##   CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## 1                 21  Provincial               Público               38
## 2                 21  Provincial               Público               38
## 3                 23     Privado               Privado               26
## 4                 23     Privado               Privado               26
## 5                 23     Privado               Privado               26
## 6                 23     Privado               Privado               26
##   PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## 1     Jujuy                  84      Susques      38084030000
## 2     Jujuy                  63    San Pedro      38063090000
## 3    Chubut                  21    Escalante      26021030018
## 4    Chubut                  35    Futaleufú      26035030000
## 5    Chubut                  35    Futaleufú      26035030000
## 6    Chubut                  35    Futaleufú      26035030000
##            LOCALIDAD   CP                             DOMICILIO          TE1
## 1            EL TORO 4411 AVENIDA LAVALLE ESQUINA 24 DE OCTUBRE             
## 2        LA MENDIETA 4522 CARLOS GARDEL S/N° - LOTE BARRO NEGRO             
## 3 COMODORO RIVADAVIA 3240                 Avenida Rivadavia 633 0297-4475310
## 4             ESQUEL 9200                        25 de Mayo 498   2945452050
## 5             ESQUEL 9200                         Ameghino 1033 02945-452155
## 6             ESQUEL 9200                          Belgrano 347             
##            TE2 TE3 TE4   ESTADO_GEO                    geometry
## 1                      SINCONFIRMAR POINT (-66.82709 -23.20853)
## 2                      SINCONFIRMAR POINT (-64.93175 -24.30814)
## 3                        CONFIRMADO POINT (-67.48355 -45.86213)
## 4 294515503809           CONFIRMADO POINT (-71.31971 -42.91505)
## 5                        CONFIRMADO POINT (-71.31803 -42.91577)
## 6                      SINCONFIRMAR POINT (-71.31992 -42.91788)

Como ven el objeto salud pasó de ser un dataframe a un sf dataframe, o sea un dataframe geográfico. Y también podemos ver todas las características del objeto. El tipo de geometría es punto, tiene 2 dimensiones y el CRS es 4326 el mismo que el de radios.

Ahora podemos mapear. Y lo hacemos igual que con el dataset de barrios y radios.

ggplot()+
  geom_sf(data=salud)

No tenemos los límites de Argentina pero claramente los puntos de los efectores parecen distribuirse en todo el país. Pero nosotros queremos trabajar solo con datos de la Ciudad de Buenos Aires! Una opción es filtrar a través de un atributo “no-geográfico” como por ejemplo el nombre de la Provincia, pero en este caso no tenemos esa información!

A no desesperarse! Una de las bondades de los datos geográficos es que podemos hacer consultas espaciales. En este caso vamos a pedirle a R que espacialmente nos filtre y seleccione todos los efectores de salud que caen dentro de algún barrio de Buenos Aires. Para eso usamos st_intersection una de las funciones espaciales de la librería sf.

Cargamos un nuevo dataframe con el polígono de la Ciudad de Buenos Aires.

caba <- read_sf("https://raw.githubusercontent.com/martoalalu/clase-geo-salud/clase-2021/data/caba.geojson")

ggplot(caba) +
  geom_sf()

Queremos quedarnos solo con los efectores de salud (puntos) que caen dentro de la Ciudad de Buenos Aires (polígono). Entonces indicamos nuestro objeto geográfico a filtrar (efectores de salud) y luego el dataframe con el que filtrar (poligono de Buenos Aires).

#Nos quedamos solo con los efectores de salud que están dentro de barrios
salud_caba <- st_intersection(salud, caba)

La consulta espacial dio como resultado un nuevo dataframe que nos indica que en la Ciudad de Buenos Aires hay 1165 efectores de salud.

Ahora sí, volvamos al mapa. Al igual que en cualquier gráfico de ggplot() podemos combinar capas. Vamos a agregar una capa el polígono de la Ciudad de Buenos Aires detrás para tener la referencia.

ggplot()+
  geom_sf(data=caba)+
  geom_sf(data=salud_caba)

A priori la distribución es parecida a la de densidad poblacional!

Ejercicio. Supongamos que querramos hacer un mapa que muestre la capacidad de respuesta del sistema de salud porteño, o sea sólo aquellos efectores que tengan internación. ¿Cómo lo hacemos? Va tip, el campo TIPOLOGÍA nos va a ayudar.

Volvamos al mapa, pero profundizando el análisis. Vamos a pintar cada efector según la columna ORIGEN_FINANCIAMIENTO el cual indica si es público o privado para ver cómo es la distribución espacial.

ggplot()+
  geom_sf(data=caba)+
  geom_sf(data=salud_caba, aes(color=ORIGEN_FINANCIAMIENTO))

Para emprolijar podemos hacer facetar el mapa, al igual que con cualquier visualización de ggplot().

ggplot()+
  geom_sf(data=caba)+
  geom_sf(data=salud_caba, aes(color=ORIGEN_FINANCIAMIENTO))+
  facet_wrap(~ORIGEN_FINANCIAMIENTO)

4 Más allá del choropleth

Hagamos otros tipos de mapas. Tener ubicaciones exactas, es decir objetos que son puntos, nos permiten hacer otros tipos de mapas como de calor (heatmap), de burbujas (bubblemap) o de densidad (densitymap), entre otros.

Dejemos de usar un rato la librería sf y vamos a considerar nuestro dataframe de efectores de salud en Buenos Aires como cualquier otro no geográfico.

Antes con st_as_coord() habíamos transformado las columnas de latitud y longitud en una geometry ahora vamos a reestablecer esas columnas, o sea vamos a hacer que dejen de ser columnas geográficas y ser sólo “numéricas”.

#Reestablecemos las columnas lat y long con sus coordenadas
salud_caba$long <- st_coordinates(salud_caba)[,1]
salud_caba$lat <- st_coordinates(salud_caba)[,2]

#Le decimos a R que salud_caba_intern deja de ser un objeto geográfico
st_geometry(salud_caba) <- NULL

class(salud_caba)

## [1] "data.frame"

Ahora salud_caba es un dataframe normal que tiene una columna latitud y otra longitud que expresan la posición de los elementos en la superficie terrestre. O sea que si hacemos un

Volvamos a ggplot() y hagamos un geom_point() tradicional.

ggplot() +
  geom_sf(data=caba) +
  geom_point(data = salud_caba, aes(x = long, y = lat,color = ORIGEN_FINANCIAMIENTO))

Ven que queda igual que si fuera un objeto geográfico! Aprovechemos nuevas funciones de ggplot() para hacer nuevos gráficos que “simulen” ser un mapa. Un mapa de densidad de puntos que muestre las zonas con mayor concentración de efectores.

ggplot() +
    geom_sf(data=caba) +
    geom_bin2d(data = salud_caba, aes(x = long, y = lat), bins = 50) +
    scale_fill_viridis_c()+
    labs(title = "Concentración de efectores de salud",
         fill = "Cantidad")

Si no nos gusta que sean cuadrados pueden ser hexagonos!

ggplot() +
    geom_sf(data=caba) +
    geom_hex(data = salud_caba, aes(x = long, y = lat), bins = 50) +
    scale_fill_viridis_c()+
    labs(title = "Concentración de efectores de salud",
         fill = "Cantidad")

O bien, las zonas de calor pueden ser más estilo “curvas de nivel”.

ggplot() +
  geom_sf(data=caba) +
  stat_density2d(data = salud_caba, aes(x = long, y = lat, fill=stat(level)),geom="polygon")+
  scale_fill_viridis_c()

5 Calculando distancias

Calcular distancias es una de las funciones espaciales más útiles ya que permite ver con claridad cómo la accesibilidad a distintos bienes y servicios se distribuye desigualmente en el espacio. En tiempos de coronavirus y cuarentena ésta cuestión cobró mayor relevancia ya que al limitarse la movilidad de las personas se cristalizan los patrones espaciales, haciendo que una parte de la población tenga un acceso privilegiado a bienes y servicios mientras que otra queda se ve obligada a trasladarse más para conseguir lo mismo (y por ende se expone más al virus).

Una de las cuestiones que surgió fue la limitación para moverse a no más de 500 metros a la redonda. Analicemos esta medida en términos de acceso a un espacio verde, algo que la Organización Mundial de la Salud considera vital para el bienestar de la población. ¿Qué porción de la población de la Ciudad de Buenos Aires está a 500 metros de una plaza?

Hagamos un choropleth map que muestre con mayor intensidad aquellos radios censales que están más cerca de una plaza.

Para esto usaremos 2 datasets: radios censales y ubicación de espacios verdes, del portal de datos abiertos del Gobierno de la Ciudad de Buenos Aires.

5.1 Explorando los datos

#Cargamos los datos
radios <- st_read('https://github.com/martoalalu/clase-geo-salud/raw/master/data/caba_radios.geojson')

## Reading layer `caba_radios' from data source `https://github.com/martoalalu/clase-geo-salud/raw/master/data/caba_radios.geojson' using driver `GeoJSON'
## Simple feature collection with 3554 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53092 ymin: -34.70574 xmax: -58.33455 ymax: -34.528
## Geodetic CRS:  WGS 84

plazas <- st_read("https://github.com/martoalalu/clase-geo-salud/raw/master/data/espacios-verdes-catastrales.geojson")

## Reading layer `espacios-verdes-catastrales' from data source `https://github.com/martoalalu/clase-geo-salud/raw/master/data/espacios-verdes-catastrales.geojson' using driver `GeoJSON'
## Simple feature collection with 1402 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53175 ymin: -34.70481 xmax: -58.33983 ymax: -34.52654
## Geodetic CRS:  WGS 84

ggplot() +
  geom_sf(data=plazas, fill="green")

Tenemos 1402 espacios verdes. Veamos qué categorías hay en el campo TIPO_EV.

unique(plazas$TIPO_EV)

##  [1] "PLAZA"                       "PLAZOLETA"                  
##  [3] "JARDÃ\u008dN"                "PARQUE"                     
##  [5] "JARDÃ\u008dN BOTÃ\u0081NICO" "PASEO"                      
##  [7] "CANTERO CENTRAL"             "CANTEROS CENTRALES"         
##  [9] NA                            "GENERAL PAZ"                
## [11] "ROTONDA"                     "DELLEPIANE"                 
## [13] "PATIO RECREATIVO NÂº12"      "PATIO PORTEÃ‘O"             
## [15] "PATIO RECREATIVO"            "CANTEROS"                   
## [17] "CANTERO"                     "JARDINES"                   
## [19] "BOSQUE"                      "ESPACIO PÃšBLICO"           
## [21] "PLAZOLETAS"                  "RESERVA ECOLÃ“GICA"         
## [23] "PATIO DE JUEGOS"             "JARDÃ\u008dN ZOOLÃ“GICO"    
## [25] "PARQUE DEPORTIVO"            "ESPACIO VERDE"

Vemos que hay espacios verdes que son canteros, plazoletas o en autopistas. Nuestro objetivo es ver el acceso a espacios verdes en los que las personas puedan distenderse y realizar algún tipo de actividad, por eso vamos a quedarnos sólo con las plazas, parques, jardín botánico, bosque, reserva ecológica, parque deportivo y espacios verdes.

plazas <- filter(plazas,TIPO_EV %in% c("PLAZA","PARQUE","JARDÍN BOTÁNICO","BOSQUE","RESERVA ECOLÓGICA","ESPACIO VERDE","PARQUE DEPORTIVO", "ESPACIO PÚBLICO"))

ggplot() +
  geom_sf(data=plazas, fill="green")

5.2 Del polígono al centroide

Para calcular distancias necesitamos hacerlo de un punto exacto a otro. Necesitamos transformar nuestros polígonos (plazas y radios) en puntos exactos, para eso vamos a crear 2 nuevas capas que tengan el centroide de cada plaza y radio. Con st_point_on_surface nos aseguramos que el punto caiga adentro del polígono.

plazas_c <- st_point_on_surface(plazas)
radios_c <- st_point_on_surface(radios)

Volvamos al mapa.

ggplot() +
  geom_sf(data=radios_c, size=0.5) +
  geom_sf(data=radios, fill=NA) +
  geom_sf(data=plazas, fill=NA) +
  geom_sf(data=plazas_c, color="lightgreen")

5.3 De un centroide a otro

Perfecto. Ahora bien, lo que queremos hacer es para cada centroide del radio censal encontrar el centroide de la plaza más cercana, calcular su distancia y agregarla como columna al dataframe de radios.

Se hace en pocas líneas pero vamos a explicarlo paso a paso

Primero tenemos que pedirle a R que por cada centroide de radio censal busque el centroide de la plaza más cercana. Esto lo hacemos utilizando la función st_nearest_feature.

st_nearest_feature(radios_c,plazas_c)[1:5]

## [1] 283  12   2 153 333

radios_c[1,]

## Simple feature collection with 1 feature and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.36622 ymin: -34.58877 xmax: -58.36622 ymax: -34.58877
## Geodetic CRS:  WGS 84
##   RADIO_ID BARRIO COMUNA POBLACION VIVIENDAS HOGARES HOGARES_NBI AREA_KM2
## 1    1_1_1 RETIRO      1       336        82      65          19 1.798997
##                      geometry
## 1 POINT (-58.36622 -34.58877)

plazas_c[283,]

## Simple feature collection with 1 feature and 26 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.37197 ymin: -34.59004 xmax: -58.37197 ymax: -34.59004
## Geodetic CRS:  WGS 84
##     OBJECTID SECCION MANZANA PARCELA         SMP TIPO_EV    NOMBRE_EV
## 283       17     003     068     000 003-068-000   PLAZA CANADÃ\u0081
##                                                                                                    UBICACION
## 283 ANTARTIDA ARGENTINA, AV. - SAN MARTIN - MARTINEZ ZUVIRIA, GUSTAVO  DR. -  RAMOS MEJIA, JOSE M., DR., AV.
##                          OBS BARRIO COMUNA SUPERFICIE LEY FECHA_LEY ORDENANZA
## 283 DENOMINACIÃ“N VERIFICADA RETIRO      1      15123   0      <NA>     16976
##      FECHA_ORD DECRETO  FECHA_DEC BOLETIN_OF   FECHA_BO FUENTE1      FUENTE2
## 283 30/12/1960   22257 23/12/1960      11505 1961-01-09    <NA> PLANO INDICE
##     FUENTE3 FUENTE4 DOMINIO NIVEL                    geometry
## 283    <NA>    <NA>    GCBA     0 POINT (-58.37197 -34.59004)

Lo que nos devuelve es una lista en la que asocia el index de cada elemento de radios_c con el index más cercano de plazas_c. Así sabemos que el primer radio censal (index = 1, Radio_id = 1_1_1) tiene como plaza más cercana a la ubicada en el index = 283 (Plaza Canadá).

Lo que vamos a hacer es calcular distancia elemento por elemento y para ello necesitamos tener “parejas” de radios y plazas, que por cada radio haya sí o sí una plaza asociada. Esto es sumamente útil ya que R va a agarrar el primer elemento de radios_c y va a calcular la distancia a su plaza asociada por orden, es decir uno a uno.

Entonces el siguiente paso es tener un dataframe que no sólo tenga el index de la plaza más cercana a cada radio sino tener la información geográfica (especialmente la columna geometry) para poder luego calcular la distancia propiamente dicha.

Vamos a traer la información de las plazas asociada al resultado que tuvimos con st_nearest_feature().

plazas_c[st_nearest_feature(radios_c, plazas_c),][1:5]

## Simple feature collection with 3554 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.52777 ymin: -34.69755 xmax: -58.35335 ymax: -34.54134
## Geodetic CRS:  WGS 84
## First 10 features:
##       OBJECTID SECCION        MANZANA    PARCELA                    SMP
## 283         17     003            068        000            003-068-000
## 12          49     005            022       024a           005-022-024a
## 2            5     001            001 018a, 020b     001-001-018a, 020b
## 153        533     001           059A  000A,000B     001-059A-000A,000B
## 333         42     005 075A,075B,075C        000 005-075A,075B,075C-000
## 333.1       42     005 075A,075B,075C        000 005-075A,075B,075C-000
## 333.2       42     005 075A,075B,075C        000 005-075A,075B,075C-000
## 12.1        49     005            022       024a           005-022-024a
## 333.3       42     005 075A,075B,075C        000 005-075A,075B,075C-000
## 206       1045     012           008C        000           012-008C-000
##                          geometry
## 283   POINT (-58.37197 -34.59004)
## 12    POINT (-58.38813 -34.60519)
## 2     POINT (-58.37828 -34.60782)
## 153   POINT (-58.36943 -34.60478)
## 333   POINT (-58.38162 -34.60371)
## 333.1 POINT (-58.38162 -34.60371)
## 333.2 POINT (-58.38162 -34.60371)
## 12.1  POINT (-58.38813 -34.60519)
## 333.3 POINT (-58.38162 -34.60371)
## 206   POINT (-58.38767 -34.60916)

nrow(plazas_c[st_nearest_feature(radios_c, plazas_c),])

## [1] 3554

Presten atención a la longitud de filas que devuelve este filtro, es 3554, la misma cantidad de observaciones que radios_c. Si bien estamos usando el dataframe plazas_c lo que le pedimos a R es que traiga los valores asociados sobre la función st_nearest_feature, que devolvía para cada index de radios_c el index de plaza_c más cercana. Asi el index = 1 del df de radios tiene como plaza más cercana al primer valor del calculo reciente (el index 283 de plazas, el object_id = 17, la plaza Canadá).

Perfecto, ya tenemos un nuevo dataframe de plazas cercanas a los radios el cual está ordenado de modo tal que el primer elemento de radios_c tiene como plaza más cercana al resultado de este filtro.

Lo que nos falta es saber la distancia en sí, lo cual hacemos con st_distance. El primer elemento de la función es el punto de origen, en este caso el centroide del radio censal, y el segundo elemento es el centroide de la plaza más cercana, el cual obtuvimos en los pasos anteriores. Luego le decimos a R que aplique la función elemento por elemento, es decir que calcule la distancia del primer elemento de radios_c al primer elemento de plazas_c[st_nearest_feature(radios_c, plazas_c),], que como sabemos es su plaza más cercana.

Veamos lo que devuelve esta función para 2 registros nada más.

st_distance(radios_c, plazas_c[st_nearest_feature(radios_c, plazas_c),], by_element = TRUE)[1:2]

## Units: [m]
## [1] 545.6843 270.0224

Perfecto! Ya tenemos el listado con las distancias calculadas desde el centroide de los radios al centroide de las plazas! Recordemos que el objetivo es hacer un choropleth map con la distancia de los radios a las plazas, y para eso necesitamos polígonos, no puntos. Necesitamos que los centroides de los radios vuelvan a ser polígonos.

Aca vamos a hacer un truquito. Como el largo del df de poligonos de radios es el mismo que el de los centroides (porque fue una transformación del primero) vamos a calcular las distancias entre el centroide de los radios y el de las plazas pero el resultado de esta operación lo vamos a agregar como columna en el dataframe de poligonos de radios censales. Como tiene el mismo largo R ni se va a enterar y nos va a dejar agregarlo sin problema!

(Acá tengan un toque de paciencia!)

radios <- radios %>% #Tomamos de base el df poligonos de radios censales
  mutate(distancia=st_distance(radios_c, plazas_c[st_nearest_feature(radios_c, plazas_c),], by_element = TRUE)) %>% #Agregamos columna nueva con los valores del calculo de distancia entre centroidoes ;)
  mutate(distancia=as.numeric(distancia)) #Convertimos a número

Miremos como quedó el dataframe de radios.

head(radios)

## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.38593 ymin: -34.60846 xmax: -58.35054 ymax: -34.57865
## Geodetic CRS:  WGS 84
##   RADIO_ID      BARRIO COMUNA POBLACION VIVIENDAS HOGARES HOGARES_NBI
## 1    1_1_1      RETIRO      1       336        82      65          19
## 2   1_12_1 SAN NICOLAS      1       341       365     116          25
## 3  1_12_10 SAN NICOLAS      1       296       629     101           1
## 4  1_12_11 SAN NICOLAS      1       528       375     136           7
## 5   1_12_2 SAN NICOLAS      1       229       445     129          16
## 6   1_12_3 SAN NICOLAS      1       723       744     314         104
##     AREA_KM2                       geometry distancia
## 1 1.79899705 MULTIPOLYGON (((-58.37189 -...  545.6843
## 2 0.01856469 MULTIPOLYGON (((-58.38593 -...  270.0224
## 3 0.04438025 MULTIPOLYGON (((-58.37879 -...  259.1591
## 4 0.36634000 MULTIPOLYGON (((-58.36733 -...  239.2970
## 5 0.01836301 MULTIPOLYGON (((-58.38454 -...  242.6532
## 6 0.03672540 MULTIPOLYGON (((-58.38154 -...  186.1420

Excelente, tenemos una columna nueva (distancia) con un número que nos va a permitir hacer nuestro querido choropleth.

Al mapa directo entonces.

ggplot() + 
  geom_sf(data = radios, aes(fill = distancia), color = NA) +
  scale_fill_viridis_c() +
  labs(title = "Distancia a plaza más cercana",
       subtitle = "Ciudad Autónoma de Buenos Aires",
       fill = "Distancia a plaza más cercana")+
  theme_void()

Genial! Si queremos podemos quedarnos sólo con los radios censales que están a menos de 500 metros.

ggplot() + 
  geom_sf(data = filter(radios,distancia < 501), aes(fill = distancia), color = NA) +
  scale_fill_viridis_c() +
  labs(title = "Radios censales a menos de 500 metros de una plaza",
       subtitle = "Ciudad Autónoma de Buenos Aires",
       fill = "Distancia a plaza más cercana")+
  theme_void()

Ok. Muy lindo el mapa pero ¿de cuántas personas estamos hablando?

radios %>% 
  st_drop_geometry() %>% 
  group_by(distancia > 500) %>% 
  summarise(total = sum(POBLACION),
            pct= total / sum(radios_c$POBLACION))

## # A tibble: 2 x 3
##   `distancia > 500`   total   pct
##   <lgl>               <dbl> <dbl>
## 1 FALSE             1849086 0.640
## 2 TRUE              1041065 0.360

El 64% de la población de Buenos Aires está a menos de 500 metros de una plaza, mientras que el 36% no.

Por último podemos hacer un histograma para ver cómo se distribuye la distancia.

ggplot(radios) +
  geom_histogram(aes(distancia))

Y lo podemos facetar por comuna…

ggplot(radios) +
  geom_histogram(aes(distancia)) +
  facet_wrap(~COMUNA, scales="free_y")

6 Yapa I: Contar puntos en polígonos

Otra operación espacial bastante común es contar puntos en poligonos. ¿Cuántas plazas hay en cada radio censal?

Usamos st_join, que es un join espacial, es decir que nos permite unir 2 dataframes a través de compartir algo en común, su ubicación. La función tiene 3 parámetros, los objetos geográficos que queremos unir y el método de join, en este caso queremos saber en qué radio están contenidas cada una de las 357 plazas entonces usamos st_within.

En primer lugar vamos a ver a qué radio censal pertenece cada plaza.

st_join(plazas_c, radios, join = st_within)

## Simple feature collection with 357 features and 35 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.52777 ymin: -34.70057 xmax: -58.35335 ymax: -34.52901
## Geodetic CRS:  WGS 84
## First 10 features:
##    OBJECTID SECCION MANZANA    PARCELA                SMP TIPO_EV
## 1         1     001     062        000        001-062-000   PLAZA
## 2         5     001     001 018a, 020b 001-001-018a, 020b   PLAZA
## 3         8     002     044        000        002-044-000   PLAZA
## 4        11     002    071B        000       002-071B-000   PLAZA
## 5        15     003     065        000        003-065-000   PLAZA
## 6        16     003    067B        000       003-067B-000   PLAZA
## 7        19     003     007        000        003-007-000   PLAZA
## 8        20     003     002        000        003-002-000   PLAZA
## 9        29     004     008  001c, 021  004-008-001c, 021   PLAZA
## 10       31     004    028B        000       004-028B-000   PLAZA
##                          NOMBRE_EV
## 1                             ROMA
## 2                     ROBERTO ARLT
## 3                          DE MAYO
## 4   PRESIDENTE JUAN DOMINGO PERÃ“N
## 5          FUERZA AÃ‰REA ARGENTINA
## 6  SALVADOR MARÃ\u008dA DEL CARRIL
## 7         DOCTOR CARLOS PELLEGRINI
## 8                         LIBERTAD
## 9           ROSARIO VERA PEÃ‘ALOZA
## 10                 CORONEL DORREGO
##                                                                               UBICACION
## 1                                AV. LEANDRO N. ALEM, LAVALLE, BOUCHARD Y TUCUMÃ\u0081N
## 2                          PREDIO EN QUE FUNCIONABA LA DIRECCION DE ASISTENCIA PÃšBLICA
## 3                                           AV. RIVADAVIA, BALCARCE, YRIGOYEN Y BOLIVAR
## 4                                     MORENO, AZOPARDO, AV. BELGRANO Y AV. PASEO COLÃ“N
## 5  DEL LIBERTADOR, AV. - SAN MARTIN - RAMOS MEJIA, JOSE M., DR., AV. - GILARDI, GILARDO
## 6           RAMOS MEJIA, JOSE M., DR., AV.- MARTINEZ ZUVIRIA, GUSTAVO, DR. - SAN MARTIN
## 7                                                       ALVEAR, AV. - ARROYO - LIBERTAD
## 8                                 LIBERTAD - ALVEAR, MARCELO T. DE - CERRITO - PARAGUAY
## 9                                                  SAN JUAN, AV. E/ PIEDRAS - CHACABUCO
## 10                                 DEFENSA - HUMBERTO I - AIETA, ANSELMO, DON - BETHLEM
##                                             OBS          BARRIO.x COMUNA.x
## 1                      DENOMINACIÃ“N VERIFICADA SAN NICOLÃ\u0081S        1
## 2  DENOMINACIÃ“N VERIFICADA / EN ACTUALIZACIÃ“N SAN NICOLÃ\u0081S        1
## 3                      DENOMINACIÃ“N VERIFICADA        MONTSERRAT        1
## 4                      DENOMINACIÃ“N VERIFICADA        MONTSERRAT        1
## 5                      DENOMINACIÃ“N VERIFICADA            RETIRO        1
## 6                      DENOMINACIÃ“N VERIFICADA            RETIRO        1
## 7                      DENOMINACIÃ“N VERIFICADA            RETIRO        1
## 8                      DENOMINACIÃ“N VERIFICADA            RETIRO        1
## 9                      DENOMINACIÃ“N VERIFICADA         SAN TELMO        1
## 10                     DENOMINACIÃ“N VERIFICADA         SAN TELMO        1
##    SUPERFICIE  LEY  FECHA_LEY                 ORDENANZA  FECHA_ORD    DECRETO
## 1   11613.283    0       <NA>                     16561 08/09/1960      19867
## 2    5699.928    0       <NA>                     25710 07/06/1971       <NA>
## 3   18447.494    0       <NA>                      <NA>       <NA> 25/05/1811
## 4   10166.412 5485 03/12/2015                      <NA>       <NA>       <NA>
## 5   26505.363    0       <NA>                     37802 24/05/1982       <NA>
## 6    4416.776    0       <NA>                      1434 21/12/1925       <NA>
## 7    1705.125    0       <NA>                      <NA>       <NA>       <NA>
## 8   11379.953    0       <NA> SE SUPONE ANTERIOR A 1836       <NA>       <NA>
## 9    7335.837    0       <NA>                     35660 02/04/1980       <NA>
## 10   2078.228    0       <NA>                      <NA> 15/05/1900       <NA>
##     FECHA_DEC BOLETIN_OF   FECHA_BO FUENTE1                         FUENTE2
## 1  17/11/1960      11479 1960-11-30    <NA>                    PLANO INDICE
## 2        <NA>      14071 1971-06-18    <NA>                    PLANO INDICE
## 3        <NA>       <NA>       <NA>    <NA> PLANO INDICE Y FICHA PARCELARIA
## 4        <NA>       4801 2016-01-15  LA LEY                            <NA>
## 5        <NA>      16790 1982-06-02    <NA>                    PLANO INDICE
## 6        <NA>       ACTA 1925-12-30    <NA>                    PLANO INDICE
## 7        <NA>       ACTA 1913-09-16    <NA>          BIBLIOTECA DE CATASTRO
## 8        <NA>       <NA>       <NA>    <NA>                FICHA PARCELARIA
## 9        <NA>      16252 1980-04-10    <NA>                    PLANO INDICE
## 10       <NA>       ACTA 1900-05-15    <NA>   BIBLIOTECA Y FICHA PARCELARIA
##    FUENTE3     FUENTE4 DOMINIO NIVEL RADIO_ID    BARRIO.y COMUNA.y POBLACION
## 1     <NA>        <NA>    GCBA     0   1_9_14 SAN NICOLAS        1       253
## 2     <NA>        <NA>    GCBA     0   1_12_9 SAN NICOLAS        1      1356
## 3     <NA>        <NA>    GCBA     0   1_14_2   MONSERRAT        1      1096
## 4     <NA>        <NA>    GCBA     0   1_14_2   MONSERRAT        1      1096
## 5     <NA>        <NA>    GCBA     0    1_7_2      RETIRO        1       770
## 6     <NA>        <NA>    GCBA     0    1_1_1      RETIRO        1       336
## 7     <NA>        <NA>    GCBA     0    1_4_6      RETIRO        1      1271
## 8     <NA>        <NA>    GCBA     0    1_8_3      RETIRO        1      1104
## 9     <NA> STREET VIEW    GCBA     0   1_26_4   SAN TELMO        1       553
## 10    <NA>        <NA>    GCBA     0   1_25_2   SAN TELMO        1       388
##    VIVIENDAS HOGARES HOGARES_NBI   AREA_KM2 distancia
## 1        435     148          26 0.08551135 113.89687
## 2        895     342          57 0.21574272 167.10782
## 3        382     226          74 0.27691508 111.75510
## 4        382     226          74 0.27691508 111.75510
## 5        281     117           0 0.14200951 135.21922
## 6         82      65          19 1.79899705 545.68428
## 7        517     486          15 0.02164026 140.50435
## 8        813     448           3 0.10244670 123.43052
## 9        291     244           0 0.01833945  78.25782
## 10       297     179           2 0.03383446  76.09029
##                       geometry
## 1  POINT (-58.36978 -34.60123)
## 2  POINT (-58.37828 -34.60782)
## 3   POINT (-58.37217 -34.6085)
## 4  POINT (-58.36863 -34.61187)
## 5   POINT (-58.37381 -34.5923)
## 6  POINT (-58.37269 -34.59084)
## 7  POINT (-58.38395 -34.59133)
## 8  POINT (-58.38342 -34.59743)
## 9  POINT (-58.37629 -34.62248)
## 10 POINT (-58.37177 -34.62048)

Esta función nos devuelve un objeto geográfico (simple feature) con el radio censal al que pertenece cada una de las plazas. Por eso tenemos 357 registros, uno por cada plaza.

Como uno de los campos es RADIO_ID podemos usar group_by y summarise para contar la cantidad de registros (plazas) que hay por cada radio censal!

radios_plazas <- plazas_c %>% 
            st_join(radios, join = st_within) %>% 
            group_by(RADIO_ID) %>% 
            count() %>% 
          st_set_geometry(NULL) # removemos la geometría porque no la vamos a usar

head(arrange(radios_plazas,desc(n)))

## # A tibble: 6 x 2
##   RADIO_ID     n
##   <chr>    <int>
## 1 14_10_1     19
## 2 2_2_3       14
## 3 8_10_8      14
## 4 13_11_1      9
## 5 13_3_2       7
## 6 14_10_2      7

El radio que más plazas tiene es el 14_10_1 y tiene 19 plazas!

Como lo que queremos hacer es un choropleth de plazas en radios ahora solo nos resta agregar la columna de cantidad de plazas al dataframe que tiene el polígono de los radios censales. Hacemos un simple left_join por radio_id y listo!

radios <- left_join(radios, radios_plazas)

Ahora sí vamos al mapa para ver la distribución espacial!

ggplot() +
  geom_sf(data=radios, aes(fill=n), color = NA) +
  scale_fill_viridis_c() +
  theme_void() +
  labs(title = "Cantidad de plazas por radio censal",
       subtitle = "Ciudad de Buenos Aires",
       fill = "Plazas por radio censal",
       caption = "Elaboración propia en base a datos de BA Data")

Listo!

Dejamos algunas ideas para practicar:

¿Cuál es el barrio que más plazas tiene? ¿Podemos calcular la cantidad de plazas por km2? ¿Cuál es el barrio que tiene más km2 de espacio verde por habitante?

7 Repaso

Y asi llegamos al final de la segunda clase donde tuvimos un pantallazo de las operaciones espaciales más comúnes:

Área -> st_area
Perímetro -> st_length
Filtrar por ubicación -> st_intersection
Convertir polígonos en puntos -> st_point_on_surface
Encontrar el objeto más cercano -> st_nearest_feature
Distancia -> st_distance
Join espacial -> st_join
Quitar la espacialidad a un df -> st_set_geometry(NULL)

Clase 2: Operaciones Espaciales

Analisis espacial con R

Martín Alalú

29/04/2021